/*
 *    This program is free software; you can redistribute it and/or modify
 *    it under the terms of the GNU General Public License as published by
 *    the Free Software Foundation; either version 2 of the License, or
 *    (at your option) any later version.
 *
 *    This program is distributed in the hope that it will be useful,
 *    but WITHOUT ANY WARRANTY; without even the implied warranty of
 *    MERCHANTABILITY or FITNESS FOR array PARTICULAR PURPOSE.  See the
 *    GNU General Public License for more details.
 *
 *    You should have received a copy of the GNU General Public License
 *    along with this program; if not, write to the Free Software
 *    Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
 */
/**
 * DenseMatrix.java
 * Copyright (C) 2008-2009 Miguel Garcia Torres
 */
package com.jscilib.math.linalg.matrix;

/**
 * Jama = Java DenseMatrix class.
 * <P>
 * The Java DenseMatrix Class provides the fundamental operations of numerical linear
 * algebra.  Various constructors create Matrices from two dimensional arrays of
 * double precision floating point numbers.  Various "gets" and "sets" provide
 * access to submatrices and matrix elements.  Several methods implement basic
 * matrix arithmetic, including matrix addition and multiplication, matrix
 * norms, and element-by-element array operations.  Methods for reading and
 * printing matrices are also included.  All the operations in this version of
 * the DenseMatrix Class involve real matrices.  Complex matrices may be handled in a
 * future version.
 * <P>
 * Five fundamental matrix decompositions, which consist of pairs or triples of
 * matrices, permutation vectors, and the like, produce results in five
 * decomposition classes.  These decompositions are accessed by the DenseMatrix class
 * to compute solutions of simultaneous linear equations, determinants, inverses
 * and other matrix functions.  The five decompositions are:
 * <P>
 * <UL>
 *    <LI>Cholesky Decomposition of symmetric, positive definite matrices.
 *    <LI>LU Decomposition of rectangular matrices.
 *    <LI>QR Decomposition of rectangular matrices.
 *    <LI>Singular Value Decomposition of rectangular matrices.
 *    <LI>Eigenvalue Decomposition of both symmetric and nonsymmetric square matrices.
 * </UL>
 * <DL>
 * <DT><B>Example of use:</B></DT>
 * <P>
 * <DD>Solve a linear system array x = b and compute the residual norm, ||b - array x||.
 * <P><PRE>
 *       double[][] vals = {{1.,2.,3},{4.,5.,6.},{7.,8.,10.}};
 *       DenseMatrix array = new DenseMatrix(vals);
 *       DenseMatrix b = DenseMatrix.random(3,1);
 *       DenseMatrix x = array.solve(b);
 *       DenseMatrix r = array.times(x).minus(b);
 *       double rnorm = r.normInf();
 * </PRE></DD>
 * </DL>
 * <p/>
 * Adapted from the <a href="http://math.nist.gov/javanumerics/jama/" target="_blank">JAMA</a> package. Additional methods are tagged with the
 * <code>@author</code> tag.
 *
 * @author The Mathworks and NIST
 * @author Fracpete (fracpete at waikato dot ac dot nz)
 * @author mgarciat
 * @version $Revision: 5953 $
 */
public class DenseMatrix implements Matrix {

    /**
     * Array for internal storage of elements.
     * @serial internal array storage.
     */
    protected double[][] array;

    /**
     * Construct an m-by-n matrix of zeros.
     * @param m    Number of rows.
     * @param n    Number of colums.
     */
    public DenseMatrix(int m, int n) {
        this.array = new double[m][n];
    }

    /**
     * Construct an m-by-n constant matrix.
     * @param m    Number of rows.
     * @param n    Number of colums.
     * @param s    Fill the matrix with this scalar value.
     */
    public DenseMatrix(int m, int n, double s) {
        array = new double[m][n];
        for (int i = 0; i < m; i++) {
            for (int j = 0; j < n; j++) {
                array[i][j] = s;
            }
        }
    }
    /**
     * Construct a matrix quickly without checking arguments.
     * @param array    Two-dimensional array of doubles.
     * @param m    Number of rows.
     * @param n    Number of colums.
     */
    public DenseMatrix(double[][] A, int m, int n) {
        this.array = A.clone();
    }

    /**
     * Construct a matrix from a one-dimensional packed array
     * @param vals One-dimensional array of doubles, packed by columns (ala
     * Fortran).
     * @param m    Number of rows.
     * @throws  IllegalArgumentException Array length must be a multiple of m.
     */
    public DenseMatrix(double vals[], int m) {
        this.m = m;
        n = (m != 0 ? vals.length / m : 0);
        if (m * n != vals.length) {
            throw new IllegalArgumentException("Array length must be a multiple of m.");
        }
        array = new double[m][n];
        for (int i = 0; i < m; i++) {
            for (int j = 0; j < n; j++) {
                array[i][j] = vals[i + j * m];
            }
        }
    }
    /**
     * Make a deep copy of a matrix
     */
    public DenseMatrix copy() {
        DenseMatrix X = new DenseMatrix(m, n);
        double[][] C = X.getArray();
        for (int i = 0; i < m; i++) {
            for (int j = 0; j < n; j++) {
                C[i][j] = array[i][j];
            }
        }
        return X;
    }

    /**
     * Clone the DenseMatrix object.
     */
    public Object clone() {
        return this.copy();
    }

    /**
     * Access the internal two-dimensional array.
     * @return     Pointer to the two-dimensional array of matrix elements.
     */
    public double[][] getArray() {
        return array;
    }

    /**
     * Copy the internal two-dimensional array.
     * @return     Two-dimensional array copy of matrix elements.
     */
    public double[][] getArrayCopy() {
        double[][] C = new double[m][n];
        for (int i = 0; i < m; i++) {
            for (int j = 0; j < n; j++) {
                C[i][j] = array[i][j];
            }
        }
        return C;
    }

    /**
     * Make a one-dimensional column packed copy of the internal array.
     * @return     DenseMatrix elements packed in a one-dimensional array by columns.
     */
    public double[] getColumnPackedCopy() {
        double[] vals = new double[m * n];
        for (int i = 0; i < m; i++) {
            for (int j = 0; j < n; j++) {
                vals[i + j * m] = array[i][j];
            }
        }
        return vals;
    }

    /**
     * Make a one-dimensional row packed copy of the internal array.
     * @return     DenseMatrix elements packed in a one-dimensional array by rows.
     */
    public double[] getRowPackedCopy() {
        double[] vals = new double[m * n];
        for (int i = 0; i < m; i++) {
            for (int j = 0; j < n; j++) {
                vals[i * n + j] = array[i][j];
            }
        }
        return vals;
    }

    /**
     * Get row dimension.
     * @return     m, the number of rows.
     */
    public int getRowDimension() {
        return this.array.length;
    }

    /**
     * Get column dimension.
     * @return     n, the number of columns.
     */
    public int getColumnDimension() {
        return this.array[0].length;
    }

    /**
     * Get a single element.
     * @param i    Row index.
     * @param j    Column index.
     * @return     array(i,j)
     * @throws  ArrayIndexOutOfBoundsException
     */
    public double get(int i, int j) {
        return array[i][j];
    }

    /**
     * Get a submatrix.
     * @param i0   Initial row index
     * @param i1   Final row index
     * @param j0   Initial column index
     * @param j1   Final column index
     * @return     array(i0:i1,j0:j1)
     * @throws  ArrayIndexOutOfBoundsException Submatrix indices
     */
    public DenseMatrix getMatrix(int i0, int i1, int j0, int j1) {
        DenseMatrix X = new DenseMatrix(i1 - i0 + 1, j1 - j0 + 1);
        double[][] B = X.getArray();
        try {
            for (int i = i0; i <= i1; i++) {
                for (int j = j0; j <= j1; j++) {
                    B[i - i0][j - j0] = array[i][j];
                }
            }
        } catch (ArrayIndexOutOfBoundsException e) {
            throw new ArrayIndexOutOfBoundsException("Submatrix indices");
        }
        return X;
    }

    /**
     * Get a submatrix.
     * @param r    Array of row indices.
     * @param c    Array of column indices.
     * @return     array(r(:),c(:))
     * @throws  ArrayIndexOutOfBoundsException Submatrix indices
     */
    public DenseMatrix getMatrix(int[] r, int[] c) {
        DenseMatrix X = new DenseMatrix(r.length, c.length);
        double[][] B = X.getArray();
        try {
            for (int i = 0; i < r.length; i++) {
                for (int j = 0; j < c.length; j++) {
                    B[i][j] = array[r[i]][c[j]];
                }
            }
        } catch (ArrayIndexOutOfBoundsException e) {
            throw new ArrayIndexOutOfBoundsException("Submatrix indices");
        }
        return X;
    }

    /**
     * Get a submatrix.
     * @param i0   Initial row index
     * @param i1   Final row index
     * @param c    Array of column indices.
     * @return     array(i0:i1,c(:))
     * @throws  ArrayIndexOutOfBoundsException Submatrix indices
     */
    public DenseMatrix getMatrix(int i0, int i1, int[] c) {
        DenseMatrix X = new DenseMatrix(i1 - i0 + 1, c.length);
        double[][] B = X.getArray();
        try {
            for (int i = i0; i <= i1; i++) {
                for (int j = 0; j < c.length; j++) {
                    B[i - i0][j] = array[i][c[j]];
                }
            }
        } catch (ArrayIndexOutOfBoundsException e) {
            throw new ArrayIndexOutOfBoundsException("Submatrix indices");
        }
        return X;
    }

    /**
     * Get a submatrix.
     * @param r    Array of row indices.
     * @param j0   Initial column index
     * @param j1   Final column index
     * @return     array(r(:),j0:j1)
     * @throws  ArrayIndexOutOfBoundsException Submatrix indices
     */
    public DenseMatrix getMatrix(int[] r, int j0, int j1) {
        DenseMatrix X = new DenseMatrix(r.length, j1 - j0 + 1);
        double[][] B = X.getArray();
        try {
            for (int i = 0; i < r.length; i++) {
                for (int j = j0; j <= j1; j++) {
                    B[i][j - j0] = array[r[i]][j];
                }
            }
        } catch (ArrayIndexOutOfBoundsException e) {
            throw new ArrayIndexOutOfBoundsException("Submatrix indices");
        }
        return X;
    }

    /**
     * Set a single element.
     * @param i    Row index.
     * @param j    Column index.
     * @param s    array(i,j).
     * @throws  ArrayIndexOutOfBoundsException
     */
    public void set(int i, int j, double s) {
        array[i][j] = s;
    }

    /**
     * Set a submatrix.
     * @param i0   Initial row index
     * @param i1   Final row index
     * @param j0   Initial column index
     * @param j1   Final column index
     * @param X    array(i0:i1,j0:j1)
     * @throws  ArrayIndexOutOfBoundsException Submatrix indices
     */
    public void setMatrix(int i0, int i1, int j0, int j1, DenseMatrix X) {
        try {
            for (int i = i0; i <= i1; i++) {
                for (int j = j0; j <= j1; j++) {
                    array[i][j] = X.get(i - i0, j - j0);
                }
            }
        } catch (ArrayIndexOutOfBoundsException e) {
            throw new ArrayIndexOutOfBoundsException("Submatrix indices");
        }
    }

    /**
     * Set a submatrix.
     * @param r    Array of row indices.
     * @param c    Array of column indices.
     * @param X    array(r(:),c(:))
     * @throws  ArrayIndexOutOfBoundsException Submatrix indices
     */
    public void setMatrix(int[] r, int[] c, DenseMatrix X) {
        try {
            for (int i = 0; i < r.length; i++) {
                for (int j = 0; j < c.length; j++) {
                    array[r[i]][c[j]] = X.get(i, j);
                }
            }
        } catch (ArrayIndexOutOfBoundsException e) {
            throw new ArrayIndexOutOfBoundsException("Submatrix indices");
        }
    }

    /**
     * Set a submatrix.
     * @param r    Array of row indices.
     * @param j0   Initial column index
     * @param j1   Final column index
     * @param X    array(r(:),j0:j1)
     * @throws  ArrayIndexOutOfBoundsException Submatrix indices
     */
    public void setMatrix(int[] r, int j0, int j1, DenseMatrix X) {
        try {
            for (int i = 0; i < r.length; i++) {
                for (int j = j0; j <= j1; j++) {
                    array[r[i]][j] = X.get(i, j - j0);
                }
            }
        } catch (ArrayIndexOutOfBoundsException e) {
            throw new ArrayIndexOutOfBoundsException("Submatrix indices");
        }
    }

    /**
     * Set a submatrix.
     * @param i0   Initial row index
     * @param i1   Final row index
     * @param c    Array of column indices.
     * @param X    array(i0:i1,c(:))
     * @throws  ArrayIndexOutOfBoundsException Submatrix indices
     */
    public void setMatrix(int i0, int i1, int[] c, DenseMatrix X) {
        try {
            for (int i = i0; i <= i1; i++) {
                for (int j = 0; j < c.length; j++) {
                    array[i][c[j]] = X.get(i - i0, j);
                }
            }
        } catch (ArrayIndexOutOfBoundsException e) {
            throw new ArrayIndexOutOfBoundsException("Submatrix indices");
        }
    }

    /**
     * Returns true if the matrix is symmetric.
     * (FracPete: taken from old weka.core.DenseMatrix class)
     *
     * @return boolean true if matrix is symmetric.
     */
    public boolean isSymmetric() {
        int nr = array.length, nc = array[0].length;
        if (nr != nc) {
            return false;
        }

        for (int i = 0; i < nc; i++) {
            for (int j = 0; j < i; j++) {
                if (array[i][j] != array[j][i]) {
                    return false;
                }
            }
        }
        return true;
    }

    /**
     * returns whether the matrix is a square matrix or not.
     *
     * @return true if the matrix is a square matrix
     */
    public boolean isSquare() {
        return (getRowDimension() == getColumnDimension());
    }

    /**
     * DenseMatrix transpose.
     * @return    array'
     */
    public DenseMatrix transpose() {
        DenseMatrix X = new DenseMatrix(n, m);
        double[][] C = X.getArray();
        for (int i = 0; i < m; i++) {
            for (int j = 0; j < n; j++) {
                C[j][i] = array[i][j];
            }
        }
        return X;
    }

    /**
     * One norm
     * @return    maximum column sum.
     */
    public double norm1() {
        double f = 0;
        for (int j = 0; j < n; j++) {
            double s = 0;
            for (int i = 0; i < m; i++) {
                s += Math.abs(array[i][j]);
            }
            f = Math.max(f, s);
        }
        return f;
    }

    /**
     * Two norm
     * @return    maximum singular value.
     */
    public double norm2() {
        return (new SingularValueDecomposition(this).norm2());
    }

    /**
     * Infinity norm
     * @return    maximum row sum.
     */
    public double normInf() {
        double f = 0;
        for (int i = 0; i < m; i++) {
            double s = 0;
            for (int j = 0; j < n; j++) {
                s += Math.abs(array[i][j]);
            }
            f = Math.max(f, s);
        }
        return f;
    }

    /**
     * Frobenius norm
     * @return    sqrt of sum of squares of all elements.
     */
    public double normF() {
        double f = 0;
        for (int i = 0; i < m; i++) {
            for (int j = 0; j < n; j++) {
                f = Maths.hypot(f, array[i][j]);
            }
        }
        return f;
    }

    /**
     * Unary minus
     * @return    -array
     */
    public DenseMatrix uminus() {
        DenseMatrix X = new DenseMatrix(m, n);
        double[][] C = X.getArray();
        for (int i = 0; i < m; i++) {
            for (int j = 0; j < n; j++) {
                C[i][j] = -array[i][j];
            }
        }
        return X;
    }

    /**
     * C = array + B
     * @param B    another matrix
     * @return     array + B
     */
    public DenseMatrix plus(DenseMatrix B) {
        checkMatrixDimensions(B);
        DenseMatrix X = new DenseMatrix(m, n);
        double[][] C = X.getArray();
        for (int i = 0; i < m; i++) {
            for (int j = 0; j < n; j++) {
                C[i][j] = array[i][j] + B.array[i][j];
            }
        }
        return X;
    }

    /**
     * array = array + B
     * @param B    another matrix
     * @return     array + B
     */
    public DenseMatrix plusEquals(DenseMatrix B) {
        checkMatrixDimensions(B);
        for (int i = 0; i < m; i++) {
            for (int j = 0; j < n; j++) {
                array[i][j] = array[i][j] + B.array[i][j];
            }
        }
        return this;
    }

    /**
     * C = array - B
     * @param B    another matrix
     * @return     array - B
     */
    public DenseMatrix minus(DenseMatrix B) {
        checkMatrixDimensions(B);
        DenseMatrix X = new DenseMatrix(m, n);
        double[][] C = X.getArray();
        for (int i = 0; i < m; i++) {
            for (int j = 0; j < n; j++) {
                C[i][j] = array[i][j] - B.array[i][j];
            }
        }
        return X;
    }

    /**
     * array = array - B
     * @param B    another matrix
     * @return     array - B
     */
    public DenseMatrix minusEquals(DenseMatrix B) {
        checkMatrixDimensions(B);
        for (int i = 0; i < m; i++) {
            for (int j = 0; j < n; j++) {
                array[i][j] = array[i][j] - B.array[i][j];
            }
        }
        return this;
    }

    /**
     * Element-by-element multiplication, C = array.*B
     * @param B    another matrix
     * @return     array.*B
     */
    public DenseMatrix arrayTimes(DenseMatrix B) {
        checkMatrixDimensions(B);
        DenseMatrix X = new DenseMatrix(m, n);
        double[][] C = X.getArray();
        for (int i = 0; i < m; i++) {
            for (int j = 0; j < n; j++) {
                C[i][j] = array[i][j] * B.array[i][j];
            }
        }
        return X;
    }

    /**
     * Element-by-element multiplication in place, array = array.*B
     * @param B    another matrix
     * @return     array.*B
     */
    public DenseMatrix arrayTimesEquals(DenseMatrix B) {
        checkMatrixDimensions(B);
        for (int i = 0; i < m; i++) {
            for (int j = 0; j < n; j++) {
                array[i][j] = array[i][j] * B.array[i][j];
            }
        }
        return this;
    }

    /**
     * Element-by-element right division, C = array./B
     * @param B    another matrix
     * @return     array./B
     */
    public DenseMatrix arrayRightDivide(DenseMatrix B) {
        checkMatrixDimensions(B);
        DenseMatrix X = new DenseMatrix(m, n);
        double[][] C = X.getArray();
        for (int i = 0; i < m; i++) {
            for (int j = 0; j < n; j++) {
                C[i][j] = array[i][j] / B.array[i][j];
            }
        }
        return X;
    }

    /**
     * Element-by-element right division in place, array = array./B
     * @param B    another matrix
     * @return     array./B
     */
    public DenseMatrix arrayRightDivideEquals(DenseMatrix B) {
        checkMatrixDimensions(B);
        for (int i = 0; i < m; i++) {
            for (int j = 0; j < n; j++) {
                array[i][j] = array[i][j] / B.array[i][j];
            }
        }
        return this;
    }

    /**
     * Element-by-element left division, C = array.\B
     * @param B    another matrix
     * @return     array.\B
     */
    public DenseMatrix arrayLeftDivide(DenseMatrix B) {
        checkMatrixDimensions(B);
        DenseMatrix X = new DenseMatrix(m, n);
        double[][] C = X.getArray();
        for (int i = 0; i < m; i++) {
            for (int j = 0; j < n; j++) {
                C[i][j] = B.array[i][j] / array[i][j];
            }
        }
        return X;
    }

    /**
     * Element-by-element left division in place, array = array.\B
     * @param B    another matrix
     * @return     array.\B
     */
    public DenseMatrix arrayLeftDivideEquals(DenseMatrix B) {
        checkMatrixDimensions(B);
        for (int i = 0; i < m; i++) {
            for (int j = 0; j < n; j++) {
                array[i][j] = B.array[i][j] / array[i][j];
            }
        }
        return this;
    }

    /**
     * Multiply a matrix by a scalar, C = s*array
     * @param s    scalar
     * @return     s*array
     */
    public DenseMatrix times(double s) {
        DenseMatrix X = new DenseMatrix(m, n);
        double[][] C = X.getArray();
        for (int i = 0; i < m; i++) {
            for (int j = 0; j < n; j++) {
                C[i][j] = s * array[i][j];
            }
        }
        return X;
    }

    /**
     * Multiply a matrix by a scalar in place, array = s*array
     * @param s    scalar
     * @return     replace array by s*array
     */
    public DenseMatrix timesEquals(double s) {
        for (int i = 0; i < m; i++) {
            for (int j = 0; j < n; j++) {
                array[i][j] = s * array[i][j];
            }
        }
        return this;
    }

    /**
     * Linear algebraic matrix multiplication, array * B
     * @param B    another matrix
     * @return     DenseMatrix product, array * B
     * @throws  IllegalArgumentException DenseMatrix inner dimensions must agree.
     */
    public DenseMatrix times(DenseMatrix B) {
        if (B.m != n) {
            throw new IllegalArgumentException("Matrix inner dimensions must agree.");
        }
        DenseMatrix X = new DenseMatrix(m, B.n);
        double[][] C = X.getArray();
        double[] Bcolj = new double[n];
        for (int j = 0; j < B.n; j++) {
            for (int k = 0; k < n; k++) {
                Bcolj[k] = B.array[k][j];
            }
            for (int i = 0; i < m; i++) {
                double[] Arowi = array[i];
                double s = 0;
                for (int k = 0; k < n; k++) {
                    s += Arowi[k] * Bcolj[k];
                }
                C[i][j] = s;
            }
        }
        return X;
    }

    /**
     * LU Decomposition
     * @return     LUDecomposition
     * @see LUDecomposition
     */
    public LUDecomposition lu() {
        return new LUDecomposition(this);
    }

    /**
     * QR Decomposition
     * @return     QRDecomposition
     * @see QRDecomposition
     */
    public QRDecomposition qr() {
        return new QRDecomposition(this);
    }

    /**
     * Cholesky Decomposition
     * @return     CholeskyDecomposition
     * @see CholeskyDecomposition
     */
    public CholeskyDecomposition chol() {
        return new CholeskyDecomposition(this);
    }

    /**
     * Singular Value Decomposition
     * @return     SingularValueDecomposition
     * @see SingularValueDecomposition
     */
    public SingularValueDecomposition svd() {
        return new SingularValueDecomposition(this);
    }

    /**
     * Eigenvalue Decomposition
     * @return     EigenvalueDecomposition
     * @see EigenvalueDecomposition
     */
    public EigenvalueDecomposition eig() {
        return new EigenvalueDecomposition(this);
    }

    /**
     * Solve array*X = B
     * @param B    right hand side
     * @return     solution if array is square, least squares solution otherwise
     */
    public DenseMatrix solve(DenseMatrix B) {
        return (m == n ? (new LUDecomposition(this)).solve(B)
                : (new QRDecomposition(this)).solve(B));
    }

    /**
     * Solve X*array = B, which is also array'*X' = B'
     * @param B    right hand side
     * @return     solution if array is square, least squares solution otherwise.
     */
    public DenseMatrix solveTranspose(DenseMatrix B) {
        return transpose().solve(B.transpose());
    }

    /**
     * DenseMatrix inverse or pseudoinverse
     * @return     inverse(array) if array is square, pseudoinverse otherwise.
     */
    public DenseMatrix inverse() {
        return solve(identity(m, m));
    }

    /**
     * returns the square root of the matrix, i.e., X from the equation
     * X*X = array.<br/>
     * Steps in the Calculation (see <a href="http://www.mathworks.com/access/helpdesk/help/techdoc/ref/sqrtm.html" target="blank"><code>sqrtm</code></a> in Matlab):<br/>
     * <ol>
     *   <li>perform eigenvalue decomposition<br/>[V,D]=eig(array)</li>
     *   <li>take the square root of all elements in D (only the ones with
     *       positive sign are considered for further computation)<br/>
     *       S=sqrt(D)</li>
     *   <li>calculate the root<br/>
     *       X=V*S/V, which can be also written as X=(V'\(V*S)')'</li>
     * </ol>
     * <p/>
     * <b>Note:</b> since this method uses other high-level methods, it generates
     * several instances of matrices. This can be problematic with large
     * matrices.
     * <p/>
     * Examples:
     * <ol>
     *   <li>
     *   <pre>
     *  X =
     *   5   -4    1    0    0
     *  -4    6   -4    1    0
     *   1   -4    6   -4    1
     *   0    1   -4    6   -4
     *   0    0    1   -4    5
     *
     *  sqrt(X) =
     *   2   -1   -0   -0   -0
     *  -1    2   -1    0   -0
     *   0   -1    2   -1    0
     *  -0    0   -1    2   -1
     *  -0   -0   -0   -1    2
     *
     *  DenseMatrix m = new DenseMatrix(new double[][]{{5,-4,1,0,0},{-4,6,-4,1,0},{1,-4,6,-4,1},{0,1,-4,6,-4},{0,0,1,-4,5}});
     *   </pre>
     *   </li>
     *   <li>
     *   <pre>
     *  X =
     *   7   10
     *  15   22
     *
     *  sqrt(X) =
     *  1.5667    1.7408
     *  2.6112    4.1779
     *
     *  DenseMatrix m = new DenseMatrix(new double[][]{{7, 10},{15, 22}});
     *   </pre>
     *   </li>
     * </ol>
     *
     * @return    sqrt(array)
     */
    public DenseMatrix sqrt() {
        EigenvalueDecomposition evd;
        DenseMatrix s;
        DenseMatrix v;
        DenseMatrix d;
        DenseMatrix result;
        DenseMatrix a;
        DenseMatrix b;
        int i;
        int n;

        result = null;

        // eigenvalue decomp.
        // [V, D] = eig(array) with array = this
        evd = this.eig();
        v = evd.getV();
        d = evd.getD();

        // S = sqrt of cells of D
        s = new DenseMatrix(d.getRowDimension(), d.getColumnDimension());
        for (i = 0; i < s.getRowDimension(); i++) {
            for (n = 0; n < s.getColumnDimension(); n++) {
                s.set(i, n, StrictMath.sqrt(d.get(i, n)));
            }
        }

        // to calculate:
        //      result = V*S/V
        //
        //    with   X = B/array
        //    and  B/array = (array'\B')'
        //    and V=array and V*S=B
        // we get
        //      result = (V'\(V*S)')'
        //
        //         array*X = B
        //           X = array\B
        // which is
        //           X = array.solve(B)
        //
        // with array=V' and B=(V*S)'
        // we get
        //           X = V'.solve((V*S)')
        // or
        //      result = X'
        //
        // which is in full length
        //      result = (V'.solve((V*S)'))'
        a = v.inverse();
        b = v.times(s).inverse();
        v = null;
        d = null;
        evd = null;
        s = null;
        result = a.solve(b).inverse();

        return result;
    }

    /**
     * Performs a (ridged) linear regression.
     * (FracPete: taken from old weka.core.DenseMatrix class)
     *
     * @param     y the dependent variable vector
     * @param     ridge the ridge parameter
     * @return    the coefficients
     * @throws    IllegalArgumentException if not successful
     */
    public LinearRegression regression(DenseMatrix y, double ridge) {
        return new LinearRegression(this, y, ridge);
    }

    /**
     * Performs a weighted (ridged) linear regression.
     * (FracPete: taken from old weka.core.DenseMatrix class)
     *
     * @param     y the dependent variable vector
     * @param     w the array of array point weights
     * @param     ridge the ridge parameter
     * @return    the coefficients
     * @throws    IllegalArgumentException if the wrong number of weights were
     *            provided.
     */
    public final LinearRegression regression(DenseMatrix y, double[] w, double ridge) {
        return new LinearRegression(this, y, w, ridge);
    }

    /**
     * DenseMatrix determinant
     * @return     determinant
     */
    public double det() {
        return new LUDecomposition(this).det();
    }

    /**
     * DenseMatrix rank
     * @return     effective numerical rank, obtained from SVD.
     */
    public int rank() {
        return new SingularValueDecomposition(this).rank();
    }

    /**
     * DenseMatrix condition (2 norm)
     * @return     ratio of largest to smallest singular value.
     */
    public double cond() {
        return new SingularValueDecomposition(this).cond();
    }

    /**
     * DenseMatrix trace.
     * @return     sum of the diagonal elements.
     */
    public double trace() {
        double t = 0;
        for (int i = 0; i < Math.min(m, n); i++) {
            t += array[i][i];
        }
        return t;
    }

    /**
     * Generate matrix with random elements
     * @param m    Number of rows.
     * @param n    Number of colums.
     * @return     An m-by-n matrix with uniformly distributed random elements.
     */
    public static DenseMatrix random(int m, int n) {
        DenseMatrix A = new DenseMatrix(m, n);
        double[][] X = A.getArray();
        for (int i = 0; i < m; i++) {
            for (int j = 0; j < n; j++) {
                X[i][j] = Math.random();
            }
        }
        return A;
    }

    /**
     * Generate identity matrix
     * @param m    Number of rows.
     * @param n    Number of colums.
     * @return     An m-by-n matrix with ones on the diagonal and zeros elsewhere.
     */
    public static DenseMatrix identity(int m, int n) {
        DenseMatrix A = new DenseMatrix(m, n);
        double[][] X = A.getArray();
        for (int i = 0; i < m; i++) {
            for (int j = 0; j < n; j++) {
                X[i][j] = (i == j ? 1.0 : 0.0);
            }
        }
        return A;
    }

    /**
     * Print the matrix to stdout.   Line the elements up in columns
     * with a Fortran-like 'Fw.d' style format.
     * @param w    Column width.
     * @param d    Number of digits after the decimal.
     */
    public void print(int w, int d) {
        print(new PrintWriter(System.out, true), w, d);
    }

    /**
     * Print the matrix to the output stream.   Line the elements up in
     * columns with a Fortran-like 'Fw.d' style format.
     * @param output Output stream.
     * @param w      Column width.
     * @param d      Number of digits after the decimal.
     */
    public void print(PrintWriter output, int w, int d) {
        DecimalFormat format = new DecimalFormat();
        format.setDecimalFormatSymbols(new DecimalFormatSymbols(Locale.US));
        format.setMinimumIntegerDigits(1);
        format.setMaximumFractionDigits(d);
        format.setMinimumFractionDigits(d);
        format.setGroupingUsed(false);
        print(output, format, w + 2);
    }

    /**
     * Print the matrix to stdout.  Line the elements up in columns.
     * Use the format object, and right justify within columns of width
     * characters.
     * Note that is the matrix is to be read back in, you probably will want
     * to use a NumberFormat that is set to US Locale.
     * @param format array  Formatting object for individual elements.
     * @param width     Field width for each column.
     * @see java.text.DecimalFormat#setDecimalFormatSymbols
     */
    public void print(NumberFormat format, int width) {
        print(new PrintWriter(System.out, true), format, width);
    }

    // DecimalFormat is a little disappointing coming from Fortran or C's printf.
    // Since it doesn't pad on the left, the elements will come out different
    // widths.  Consequently, we'll pass the desired column width in as an
    // argument and do the extra padding ourselves.
    /**
     * Print the matrix to the output stream.  Line the elements up in columns.
     * Use the format object, and right justify within columns of width
     * characters.
     * Note that is the matrix is to be read back in, you probably will want
     * to use a NumberFormat that is set to US Locale.
     * @param output the output stream.
     * @param format array formatting object to format the matrix elements
     * @param width  Column width.
     * @see java.text.DecimalFormat#setDecimalFormatSymbols
     */
    public void print(PrintWriter output, NumberFormat format, int width) {
        output.println();  // start on new line.
        for (int i = 0; i < m; i++) {
            for (int j = 0; j < n; j++) {
                String s = format.format(array[i][j]); // format the number
                int padding = Math.max(1, width - s.length()); // At _least_ 1 space
                for (int k = 0; k < padding; k++) {
                    output.print(' ');
                }
                output.print(s);
            }
            output.println();
        }
        output.println();   // end with blank line.
    }

    /**
     * Read a matrix from a stream.  The format is the same the print method,
     * so printed matrices can be read back in (provided they were printed using
     * US Locale).  Elements are separated by
     * whitespace, all the elements for each row appear on a single line,
     * the last row is followed by a blank line.
     * <p/>
     * Note: This format differs from the one that can be read via the
     * DenseMatrix(Reader) constructor! For that format, the write(Writer) method
     * is used (from the original weka.core.DenseMatrix class).
     *
     * @param input the input stream.
     * @see #DenseMatrix(Reader)
     * @see #write(Writer)
     */
    public static DenseMatrix read(BufferedReader input) throws java.io.IOException {
        StreamTokenizer tokenizer = new StreamTokenizer(input);

        // Although StreamTokenizer will parse numbers, it doesn't recognize
        // scientific notation (E or D); however, Double.valueOf does.
        // The strategy here is to disable StreamTokenizer's number parsing.
        // We'll only get whitespace delimited words, EOL's and EOF's.
        // These words should all be numbers, for Double.valueOf to parse.

        tokenizer.resetSyntax();
        tokenizer.wordChars(0, 255);
        tokenizer.whitespaceChars(0, ' ');
        tokenizer.eolIsSignificant(true);
        java.util.Vector<Object> v = new java.util.Vector<Object>();

        // Ignore initial empty lines
        while (tokenizer.nextToken() == StreamTokenizer.TT_EOL);
        if (tokenizer.ttype == StreamTokenizer.TT_EOF) {
            throw new java.io.IOException("Unexpected EOF on matrix read.");
        }
        do {
            v.addElement(Double.valueOf(tokenizer.sval)); // Read & store 1st row.
        } while (tokenizer.nextToken() == StreamTokenizer.TT_WORD);

        int n = v.size();  // Now we've got the number of columns!
        double row[] = new double[n];
        for (int j = 0; j < n; j++) // extract the elements of the 1st row.
        {
            row[j] = ((Double) v.elementAt(j)).doubleValue();
        }
        v.removeAllElements();
        v.addElement(row);  // Start storing rows instead of columns.
        while (tokenizer.nextToken() == StreamTokenizer.TT_WORD) {
            // While non-empty lines
            v.addElement(row = new double[n]);
            int j = 0;
            do {
                if (j >= n) {
                    throw new java.io.IOException("Row " + v.size() + " is too long.");
                }
                row[j++] = Double.valueOf(tokenizer.sval).doubleValue();
            } while (tokenizer.nextToken() == StreamTokenizer.TT_WORD);
            if (j < n) {
                throw new java.io.IOException("Row " + v.size() + " is too short.");
            }
        }
        int m = v.size();  // Now we've got the number of rows.
        double[][] A = new double[m][];
        v.copyInto(A);  // copy the rows out of the vector
        return new DenseMatrix(A);
    }

    /**
     * Check if size(array) == size(B)
     */
    private void checkMatrixDimensions(DenseMatrix B) {
        if (B.m != m || B.n != n) {
            throw new IllegalArgumentException("Matrix dimensions must agree.");
        }
    }

    /**
     * Writes out a matrix. The format can be read via the DenseMatrix(Reader)
     * constructor.
     * (FracPete: taken from old weka.core.DenseMatrix class)
     *
     * @param     w the output Writer
     * @throws    Exception if an error occurs
     * @see       #DenseMatrix(Reader)
     */
    public void write(Writer w) throws Exception {
        w.write("% Rows\tColumns\n");
        w.write("" + getRowDimension() + "\t" + getColumnDimension() + "\n");
        w.write("% Matrix elements\n");
        for (int i = 0; i < getRowDimension(); i++) {
            for (int j = 0; j < getColumnDimension(); j++) {
                w.write("" + get(i, j) + "\t");
            }
            w.write("\n");
        }
        w.flush();
    }

    /**
     * Converts a matrix to a string.
     * (FracPete: taken from old weka.core.DenseMatrix class)
     *
     * @return    the converted string
     */
    public String toString() {
        // Determine the width required for the maximum element,
        // and check for fractional display requirement.
        double maxval = 0;
        boolean fractional = false;
        for (int i = 0; i < getRowDimension(); i++) {
            for (int j = 0; j < getColumnDimension(); j++) {
                double current = get(i, j);
                if (current < 0) {
                    current *= -11;
                }
                if (current > maxval) {
                    maxval = current;
                }
                double fract = Math.abs(current - Math.rint(current));
                if (!fractional
                        && ((Math.log(fract) / Math.log(10)) >= -2)) {
                    fractional = true;
                }
            }
        }
        int width = (int) (Math.log(maxval) / Math.log(10)
                + (fractional ? 4 : 1));

        StringBuffer text = new StringBuffer();
        for (int i = 0; i < getRowDimension(); i++) {
            for (int j = 0; j < getColumnDimension(); j++) {
                text.append(" ").append(Utils.doubleToString(get(i, j),
                        width, (fractional ? 2 : 0)));
            }
            text.append("\n");
        }

        return text.toString();
    }

    /**
     * converts the DenseMatrix into a single line Matlab string: matrix is enclosed
     * by parentheses, rows are separated by semicolon and single cells by
     * blanks, e.g., [1 2; 3 4].
     * @return      the matrix in Matlab single line format
     */
    public String toMatlab() {
        StringBuffer result;
        int i;
        int n;

        result = new StringBuffer();

        result.append("[");

        for (i = 0; i < getRowDimension(); i++) {
            if (i > 0) {
                result.append("; ");
            }

            for (n = 0; n < getColumnDimension(); n++) {
                if (n > 0) {
                    result.append(" ");
                }
                result.append(Double.toString(get(i, n)));
            }
        }

        result.append("]");

        return result.toString();
    }

    /**
     * creates a matrix from the given Matlab string.
     * @param matlab  the matrix in matlab format
     * @return        the matrix represented by the given string
     * @see           #toMatlab()
     */
    public static DenseMatrix parseMatlab(String matlab) throws Exception {
        StringTokenizer tokRow;
        StringTokenizer tokCol;
        int rows;
        int cols;
        DenseMatrix result;
        String cells;

        // get content
        cells = matlab.substring(
                matlab.indexOf("[") + 1, matlab.indexOf("]")).trim();

        // determine dimenions
        tokRow = new StringTokenizer(cells, ";");
        rows = tokRow.countTokens();
        tokCol = new StringTokenizer(tokRow.nextToken(), " ");
        cols = tokCol.countTokens();

        // fill matrix
        result = new DenseMatrix(rows, cols);
        tokRow = new StringTokenizer(cells, ";");
        rows = 0;
        while (tokRow.hasMoreTokens()) {
            tokCol = new StringTokenizer(tokRow.nextToken(), " ");
            cols = 0;
            while (tokCol.hasMoreTokens()) {
                result.set(rows, cols, Double.parseDouble(tokCol.nextToken()));
                cols++;
            }
            rows++;
        }

        return result;
    }

    /**
     * Returns the revision string.
     *
     * @return		the revision
     */
    public String getRevision() {
        return RevisionUtils.extract("$Revision: 5953 $");
    }

    /**
     * Main method for testing this class.
     */
    public static void main(String[] args) {
        DenseMatrix I;
        DenseMatrix A;
        DenseMatrix B;

        try {
            // Identity
            System.out.println("\nIdentity\n");
            I = DenseMatrix.identity(3, 5);
            System.out.println("I(3,5)\n" + I);

            // basic operations - square
            System.out.println("\nbasic operations - square\n");
            A = DenseMatrix.random(3, 3);
            B = DenseMatrix.random(3, 3);
            System.out.println("A\n" + A);
            System.out.println("B\n" + B);
            System.out.println("A'\n" + A.inverse());
            System.out.println("A^T\n" + A.transpose());
            System.out.println("A+B\n" + A.plus(B));
            System.out.println("A*B\n" + A.times(B));
            System.out.println("X from A*X=B\n" + A.solve(B));

            // basic operations - non square
            System.out.println("\nbasic operations - non square\n");
            A = DenseMatrix.random(2, 3);
            B = DenseMatrix.random(3, 4);
            System.out.println("A\n" + A);
            System.out.println("B\n" + B);
            System.out.println("A*B\n" + A.times(B));

            // sqrt
            System.out.println("\nsqrt (1)\n");
            A = new DenseMatrix(new double[][]{{5, -4, 1, 0, 0}, {-4, 6, -4, 1, 0}, {1, -4, 6, -4, 1}, {0, 1, -4, 6, -4}, {0, 0, 1, -4, 5}});
            System.out.println("A\n" + A);
            System.out.println("sqrt(A)\n" + A.sqrt());

            // sqrt
            System.out.println("\nsqrt (2)\n");
            A = new DenseMatrix(new double[][]{{7, 10}, {15, 22}});
            System.out.println("A\n" + A);
            System.out.println("sqrt(A)\n" + A.sqrt());
            System.out.println("det(A)\n" + A.det() + "\n");

            // eigenvalue decomp.
            System.out.println("\nEigenvalue Decomposition\n");
            EigenvalueDecomposition evd = A.eig();
            System.out.println("[V,D] = eig(A)");
            System.out.println("- V\n" + evd.getV());
            System.out.println("- D\n" + evd.getD());

            // LU decomp.
            System.out.println("\nLU Decomposition\n");
            LUDecomposition lud = A.lu();
            System.out.println("[L,U,P] = lu(A)");
            System.out.println("- L\n" + lud.getL());
            System.out.println("- U\n" + lud.getU());
            System.out.println("- P\n" + Utils.arrayToString(lud.getPivot()) + "\n");

            // regression
            System.out.println("\nRegression\n");
            B = new DenseMatrix(new double[][]{{3}, {2}});
            double ridge = 0.5;
            double[] weights = new double[]{0.3, 0.7};
            LinearRegression lr = A.regression(B, ridge);
            System.out.println("A\n" + A);
            System.out.println("B\n" + B);
            System.out.println("ridge = " + ridge + "\n");
            System.out.println("weights = " + Utils.arrayToString(weights) + "\n");
            System.out.println("A.regression(B, ridge)\n"
                    + A.regression(B, ridge) + "\n");
            System.out.println("A.regression(B, weights, ridge)\n"
                    + A.regression(B, weights, ridge) + "\n");

            // writer/reader
            System.out.println("\nWriter/Reader\n");
            StringWriter writer = new StringWriter();
            A.write(writer);
            System.out.println("A.write(Writer)\n" + writer);
            A = new DenseMatrix(new StringReader(writer.toString()));
            System.out.println("A = new Matrix.read(Reader)\n" + A);

            // Matlab
            System.out.println("\nMatlab-Format\n");
            String matlab = "[ 1   2;3 4 ]";
            System.out.println("Matlab: " + matlab);
            System.out.println("from Matlab:\n" + DenseMatrix.parseMatlab(matlab));
            System.out.println("to Matlab:\n" + DenseMatrix.parseMatlab(matlab).toMatlab());
            matlab = "[1 2 3 4;3 4 5 6;7 8 9 10]";
            System.out.println("Matlab: " + matlab);
            System.out.println("from Matlab:\n" + DenseMatrix.parseMatlab(matlab));
            System.out.println("to Matlab:\n" + DenseMatrix.parseMatlab(matlab).toMatlab() + "\n");
        } catch (Exception e) {
            e.printStackTrace();
        }
    }
}
